Inferencia Causal en R
2024-06-07
Conceptos Iniciales
Inferencia Causal
Cuando hablamos de un proceso de inferencia causal nos referimos al proceso de hacer afirmaciones sobre una \(población\) basándonos en una \(muestra\) representativa compuesta de \(unidades\) (misma que dependerá de diferentes supuestos), particularmente, usaremos estos dos conceptos para entender el efecto de un fenómeno ocurrido de la \(muestra\) sobre la \(población\), buscando si es \(estadísticamente \ significativo\) para establecer una relación \(causal\).
Modelo de Rubin, Holland
Una afirmación ya bastante conocida es que correlación no implica causalidad, sin embargo no resulta muy claro el camino a seguir estadísticamente para crear un modelo causal, una parte importante es establecer el concepto de \(aleatoreidad\) lo cual obedece a un $proceso generador de datos $ el cual puede se puede aislar bajo sus propias condiciones.
Retomando los conceptoa anteriores, consideraremos una \(población\), \(muestra\) y \(unidades\), de las últimas podremos establecer \(funcione\) que traduciremos como \(variables\).
En ese sentido, podemos definir los componentes del modelo como:
- Unidad : Un objeto \(i\) en un punto del tiempo específico.
- Tratamiento : Una intervención \(t\), cuyos efectos queremos estimar comparando con un \(control / c\) determinado por la ausencia de intervención.
- Estado del mundo: \(S\) para la unidad \(i\) la cual tiene un estado de tratamiento definido por \(S \in {t,c}\).
- Resultado potencial: Antes de que el tratamiento es asignado, podemos estimar la medición de nuestro resultado de interés \(Y\) para cada unidad \(i\), bajo ambos tratamiento y control, identificados como \(Y_{i}^{t}\) y \(Y_{i}^{c}\).
- Mecanismo de asignación: El algoritmo que determina si cada \(i\) recibe \(t\) o \(c\).
- Efecto causal: La diferencia entre el resultado potencial para una unidad dada: \(Y_{i}^{t} - Y_{i}^{c} = T_{i}\).
Una vez que \(S = c\) podemos obtener el \(efecto /del tratamiento\) como \(Y_{i}^{t} - Y_{i}^{c} = T_{i}\). No obstante una vez que \(S\) está determinado solo observamos un resultado, nunca podemos observar dos resultados potenciales para la misma unidad, a esto se refiere el Problema Fundamental de la Inferencia Causal, es decir, es imposible observar el valor de \(Y_{i}^{t}\) y \(Y_{i}^{c}\) para la misma unidad y por tanto, es imposible observar el efecto de \(t\) contra \(c\) para una \(i\) en particular.
//
Una solución al problema anterior es, en lugar de concentrarnos en calcular \(T_{i}\) calculamos \(T\) que será el \(Efecto / medio / de/ tratamientoo\) calculado como: $T = E(Y_{i}^{t} - Y_{i}^{c}) = E(T_{i}) $ el cual se obtiene directamente recordando que la media es un operador lineal.
Aunque no siempre podemos observar \(T\), podemos estimarlo de la siguiente forma:
\(\hat{T} = \hat{E}(Y_{i}^{t}| S=t) - \hat{E}( Y_{i}^{c}|S=c)\)
Generalmente \(\hat{T}\) es un buen estimador de \(T\) (no quiere decir que sea el único), lo cual depende del mecanismo de asignación.
Consideremos: \(E(Y_{i}^{t}) = (Y_{i}^{t}| S=t)\) y \(E(Y_{i}^{c}) = (Y_{i}^{c}| S=c)\)
Entonces \(E(\hat{T}) = T\) se cumplirá si \(S\) es independeinte de \(Y\), entonces el grupo tratado es una muestra aleatoria del grupo de los potencialmente tratados, mismo caso para el grupo de control.
Notemos entonces que la asignación aleatoria de \(S\) hace que la independencia sea plausible pero no testeable.
Problema de identificación.
Cuando el supuesto de independencia entre \(S\) y \(Y\) cuando hay alguna relación a priori entre ambos, por ejemplo, si el tratamiento es parte de un procceso de optimización que involucra al resultado, caso que es común.
A partir de lo anterior, nos referimos al problema de identificación a la dificultad de obtener un estimador insesgado para identificar parámetros causales a menos que el mecanismo de asignación cumpla con distintas condiciones.
Organización y estructura de datos.
Vectores
Los vectores se usan regularmente para establecer una secuencia numérica, cada elemento del vector revela información sobre los datos en general, se representan comunmente como: \((S_{1}.... S_{6})\), en este caso se tratará de un vector de 6 elementos.
Matrices
Las matrices son objectos relacionados en dos dimensiones los cuales crean combinaciones de fila-columna, se representan comunmente como:
\[\begin{equation} \begin{pmatrix} a_{a} & b_{a} & c_{a} & d_{a} \\ a_{b} & b_{b} & c_{b} & d_{b} \\ a_{c} & b_{c} & c_{c} & d_{c} \\ a_{d} & b_{d} & c_{d} & d_{d} \\ \end{pmatrix} \end{equation}\] =
\[\begin{equation} \begin{pmatrix} 1 & 2 & 3 & 4 \\ 8 & 7 & 6 & 5 \\ 9 & 10 & 11 & 12 \\ 13 & 14 & 15 & 16 \\ \end{pmatrix} \end{equation}\]
Data frame
Un data frame contendrá información sobre múltiples unidades para diferentes vectores, los elemtnos que contienen este data frame regularmente a nivel de fila tendrán identificadores, aunque no siempre puede darse el caso en el que la información esté idenxada, esta última característica será útil si queremos agrupar nuestros datos.
Tipos de variables
Variables binarias: sólo pueden tomar dos valores de un conjunto discreto.
Variables categóricas: capturan si la unidad del dato está clasificada en un rango de categorías.
Variables ordinarias: capturan información en un rango particular.
Variables continuas: capturan datos en un margen más amplio, generalmente son datos numéricos en un conjunto continuo.
Obtención de datos
Descargas directas: generalmente son datos con una estructura definida, son de fácil acceso a través de links descargables.
Scrapeo de datos: es un proceso de automatización comunmente realizado en R o Python, en el cual se obtienen datos de cada link en una página o bien se retoman usando la estructura HTML de la página en cuestión.
Data Capture: son datos que se encuentran previamente clasificados, generalmente siguen una estructura estadística.
API: es una fuente en la que podemos acceder más directamente a los datos a través de código, generalmente los datos se encuentran en foramto JSON.
Teorema del Límite Central
Muestra aleatoria: se refiere a una realización aleatoria de un conjunto de objetos en un data frame, esperamos que al tomar subconjuntos de dicha muestra:
La media de ese subconjunto converga (en media) a la media de la muestra completa.
La varianza de ese subconjunto sea representativa a la varianza de nuestra muestra completa.
Lo anterior se cumplirá sí y sólo sí los datos son independientes e idénticamente distribuidos (no hay ninguna conexión entre las observaciones o filas para un dataframe), es decir, si no hay subgrupos en la población que podrían tener resultados correlacionados.
Si quisíeramos hacer algún tipo de inferencia partiendo de los grupos de un data frame tenemos que usar el método de aleatorización estratificada para asegurarnos de considerar submuestras representativas para cada subgrupo.
Distribuciones: nos muetran el comportamiento las varibales aleatorias, y son particularmente útiles para realizar pruebas de hipótesis. Podemos distinguir entre:
Discretas: Involucran descontinuidades escalonadas en la densidad, también tienen la naturaleza de mostrar realizaciones binarias, es por eso que los ejemplos más comunes son las distribuciones Bernoulli y Binomial.
Continuas: Tienen una forma suavizada, son un caso de una distribución discreta, la forma más común es la distribución normal (Gaussiana)
Lo anterior nos sierve para establecer el Teorema del Límite Central, el cual establece que conforme el número de realizaciones aumenta, la media muestral tiene a distribuirse de forma normal alrederor la media de la población, también la desviación estandar tiende a disminuir conforme \(n\) disminuye (a causa de la convergencia en media).
Para una realización aleatoria podemos hacer inferencias sobre qué podría ser representativo de nuestra población, conforme agregamos más datos, nuestra inferencia se vuelve más robusta.
knitr::include_graphics("Sample_subset.png")
Analítica Descriptiva: Medidas de dispersión y localización
Medidas de localización/tendencia central
Consideremos que la localización central de un data frame (para un vector) requiere considerar diferentes medidas de tendencia central como:
Media: También llamada Esperanza o Valor Esperado y denotado comunmente como \(\bar{X}\) es igual a la suma del vector dividido entre su longitud (notemos que este valor indirectamente que el peso de cada valor en el vector es el mismo en términos de probabilidad).
Mediana: El valor que se encuentra justo en medio del vector listado (el número de valores previos y posteriores es el mismo aproximádamente).
Moda: El valor con más ocurrencia en el vector.
Medidas de dispersión
Generalmente los datos están distribuidos alrededor de una tendencia central, para organizarlos mejor distinguimos entre:
Rango: Diferencia entre el valor más pequeño y el más grande en una secuencia numérica (vector numérico).
Percentiles/Quartiles: Muestran al densidad de la distribución sobre un punto en particular.
Rango IQ: Distancia entre dos segmentos de datos.
Varianza: Suma de los datos estandarizados (restándoles la media) divivido entre \(n-1\), es decir: \(\sigma^2 = 1/n-1 \sum_{1}^n (X_{i} - \bar{X})^2\) recordemos que la Desviación estándar es la raíz cuadrada de la Varianza (i.e \(\sigma\)).
Correlación
- La correlación nos ayuda a evaluar la relación entre dos variables (\(X\), \(Y\); o más), la visualización más común es el gráfico de dispersión donde podemos observar como cuando un vector aumenta su valor, los valores del otro vector aumentan (relación directa), disminuyen (relación inversa) o se quedan igual (no hay relación), podemos obtener la correlación calculando \(\sum_{1}^n (X_{i} - \bar{X})(Y_{i} - \bar{y})/ SD^X SD^Y\), esta medida deriva de la covarianza y es importante mencionar que la correlación no implica causalidad.
Visualizaciones
Tablas
Boxplots
Histogramas
Densidad
Randomized Control Trials (RCTs)
Partiendo del modelo de Rubin, podemos recordar que si el tratamiento es determinado de forma aleatoria, entonces la independencia en la muestra es plausible (mas no garantizada), los RCTs hacen uso de este supuesto para generar modelos de asignación aleatoria. Algo importante a diferencias son los tratamientos exógenos vs los endógenos.
La idea detrás es que cuando los agentes eligen el estado del mundo \(S\) (endógeno) dependederá en los resultados \(Y\).
Por otro lado, si es exógeno podemos definir un “agente optimizador fuera del control” que hace la independencia plausible.
Lo anterior se puede entender también en un problema de exogeneidad sobre una regesión continua, lo cual le resta robustez y significancia estadística a nuestro modelo y por tanto a la intepretación causal de los resultados, cuando hacemos un RCT la idea principal es resolver el problema de endogenidad teniendo una muestra independiente a priori, aunque también podemos encontrarnos con el problema de una aleatorización pobre, una forma de verificar que no ocurre dicho problema es realizando balance tests, los cuales buscan verificar la independencia del estado del mundo \(S\) con las características exógenas (no el resultado \(Y\)).
Limitantes de los RCT
Beneficios:
- Dado que podemos establecer independencia a priori, el problema de identificación se resuelve fácilmente, pudiendo establecer un efecto causal.
Limitantes o problemas potenciales:
- Aleatorización fallida
- Desgaste (Attrition)
- Efectos en equilibrio general
- Practicidad
- Validación externa
- Sesgo por elección
Attrition: Algunos participantes pueden dejar el estudio, el problema está en que no podremos encontrar \(Y\) de las personas que abandonan el estudio, si el attrition no está correlacionado \(Y^{s}\), no será un problema, dado que aún se podrá tener independencia de \(S\) y \(Y^{t}\) y \(Y^{c}\) en el resto de la muestra, si el attrition está correlacionado con \(Y^{S}\) y el tratamiento cambia el attrition entonces la indepedencia se rompe, necesitaríamos:
\(E[Y^{S}|S = s, observada] = E[Y^{S}, observada]\)
Observaríamos un conjunto no aleatorio de los resultados tratados.
Compliance y la intención a ser tratado (ITT)
En algunas ocasiones también las personas se pueden negar a tomar el tratamiento, también puede darse el casp que las personas del grupo de control busquen tratarse fuera del programa, si es que se observan los resultados, se puede obtener un estimador válido del parámetro sobre la intención a ser tratado (ITT), podemos entender mejor lo anterior a través de la siguiente Matrix de Compliance:
| Assignación t | Assignación c | |
|---|---|---|
| Recibo t | A | B |
| Recibo c | C | D |
En este caso, A y D son compliers, mientras que B y C no lo serán. De inicio, podríamos pensar que nuna buena opción sería la de comparar A y B contra C y D, pero se perdiría la independencia dado que tenemos agentes endógenos que eligen en el estado del tratamiento, esto porque la asignación de \(t\) o \(c\) no necesariamente implica que hayan recibido \(t\) o \(c\) necesariamente.
De tal forma se desarrolla el concepto de la intención a ser tratado, el cual compara A y C así como B y D, si todos los agentes son compliers entonces se cumplirá que ITT = ATT.
Efectos en Equilibrio General
Ocurre cuando el grupo de control es afectado por el tratamiento, puede darse debido a efectos en equilibrio general o bien a efectos de pares (peer effects), en caso que se identifiquen esos efectos de forma potencial, el RCT debe de ser rediseñado para prevenirles, lo cual regularmente reduce el tamaño de la muestra efectiva a través de clusters.
Practicidad y validación.
Ética: Hay experimentos que pueden generar estrés financiero o una situación poco deseable para la muestra.
Factibilidad: Existe una limitante sobre los recursos que podemos usar y el escenario que podemos plantear.
Costo: Los RCT’s tienden a ser costosos en comparación con los estudios observacionales.
Validación externa: Una estimador será válido internamente si es un estimador insesgado para la población del experimento, también será válido externamente si es un estimador insesgado para la población fuera del experimento, no todos los RCT tienen validación externa pero se busca que la tengan.
Los RCT’s serán nuestro marco estándar para estimar una inferencia causal a causa de que la independencia es plausible al controlar aleatoriamente el mecanismo de asignación, resolviendo el problema de identificación, por eso cuando examinamos datos observacionales (no experimentales) siepre es buena idea comparar con el RCT inicial para comprender los problemas potenciales, aunque es importante tener en mente las limitantes de los RCTs, teniendo claro el tipo de estimación que queremos obtener, así como los parámetros que son creíbles para estimar, considerando attrition, compliance, equilibrio general, etc.
Regresión Bivariada
También llamada OLS(Ordinary Least Squares) Bivariada o Mínimos Cuadrados Ordinarios Bivariado, nos ofrece un modelo en el que se ajusta una línea entre dos variables X y Y, por lo que se asume de inició una relación lineal entre X y Y tal que:
\[\begin{equation} Y = \beta_{0} + \beta_{1}X + u \end{equation}\]
MCO (OLS) será el método para estimar la pendiente \(\beta_{1}\) y el intercepto \(\beta_{0}\), aunque también se podrá estimar a través de máxima verosimilitud.
La ecuación estándar de regresión para una población será:
\[\begin{equation} Y_{i} = \beta_{0} + \beta_{1}X_{i} + u_{i} \end{equation}\]
Donde:
- \(Y_{i}\) es la variable dependiente.
- \(X_{i}\) es la variable independiente.
- \(\beta_{0}\) es la constante/intersecto verdadera.
- \(\beta_{1}\) es la pendiente verdadera.
- \(u_{i}\) es el término de error verdadero.
Por otro lado, nuestros estimadores serán marcados con un “sombrero” de la siguiente forma:
- \(\hat{\beta}\) será el estimador del coeficiente.
- \(\hat{Y_{i}}\) será el valor predecido.
Por lo que \(\hat{Y} = \hat{\beta_{0}} + \hat{\beta_{1}} X_{i}\)
- \(\hat{u_{i}}\) es un residuo.
Por lo que $ = Y_{i} - = Y_{i} - - X_{i}$
Observamos \(X\) y \(Y\) y estimamos \(\beta\) y \(u\).
Mínimos Cuadrados Ordinarios (MCO)
Sigue una lógica de optimización en la que se minimizan las desviaciones de los valores predecidos con los valores reales en los datos, generalmente toma una muestra de tamaño \(n\) dado por:
\[\begin{equation} min_{b_{o}b_{1}} = \sum_{i}^{n}(Y_{i} - b_{0} - b_{1}X_{i})^2 \end{equation}\]
Los residuos se refieren a los errores que el modelo está haciendo, los cuales están dados por \(\bar{u_{i}} = Y_{i} - \hat{\beta_{0}} - \hat{\beta_{1}} X_{i}\), podemos reescribir el problema de MCO como:
\[\begin{equation} min_{b_{o}b_{1}} \sum_{i=1}^{n} (\tilde{u_{i}})^2 \end{equation}\]
donde \(\tilde{u_{i}} = Y_{i} - b_{0} - b_{1}X_{i}\), la solución dada por MCO minimiza la suma de los residuos al cuadrado, es decir mnimiza la distancia vertical entre los puntos.
knitr::include_graphics("OLS.png")
Notemos sin embargo que \(\hat{\beta_{0}}\) y \(\hat{\beta_{1}}\) son los valroes de \(\beta_{0}\) y \(\beta_{1}\) que minimizan la suma de residuos al cuadrado, donde cada uno está dado por:
\[\begin{equation} \hat{\beta_{1}} = \frac{\sum_{i=1}^{n}(X_{i} - \bar{X}) (Y_{i} - \bar{Y})} {{\sum}_{i=1}^{n} (X_{i} - \bar{X})^2} \end{equation}\]
\[\begin{equation} = \frac{\hat{Cov(X,Y)}}{\hat{Var(X)}} = \frac{S_{XY}}{S_{X}^2} \end{equation}\]
\[\begin{equation} \hat{\beta_{0}} = \bar{Y} - \hat{\beta_{1}}\bar{X} \end{equation}\]
Donde
\(\hat{X} = \frac{1}{n} \sum_{i=1}^{n} X_{i}\) es la media muestral de X.
y \(\hat{Y} = \frac{1}{n} \sum_{i=1}^{n} Y_{i}\) es la media muestral de Y.
El método de MCO nos ayuda a encontrar \(\hat{\beta_{s}}\), lo cual define el valor predecido y los residuos.
Los valores predecidos (\(\hat{Y}\) es el modelo estimado de \(Y_{i}\) dado \(X_{i}\)):
\[\begin{equation} \hat{Y_{i}} = \hat{\beta_{0}} + \hat{\beta_{1}}X_{i} \end{equation}\]
Predice los valores a partir de los estimadores de MCO usando medias condicionales:
\[\begin{equation} E(Y_{i} | X_{i}) = \beta_{0} + \beta_{1}X_{i} \end{equation}\]
\[\begin{equation} \hat{E}(Y_{i} | X_{i}) = \hat{\beta_{0}} + \hat{\beta_{1}}X_{i} = \hat{Y_{i}} \end{equation}\]
Los valores predecidos son las mejores aproximaciones de \(Y_{i}\) cuando sabemos \(X_{i}\).
Los residuos serán la diferencia entre los resultados veraderos y los valores predicidos.
\[\begin{equation} \hat{u_{i}} = Y_{i} - \hat{Y_{i}} = Y_{i} - \hat{\beta_{0}} - \hat{\beta_{1}}X_{i} \end{equation}\]
Supuestos importantes para MCO - Teorema de Gauss - Markov
MCO será un “buen” método para estimar una relación causal entre \(X\) y \(Y\) a través del estimador/pendiente \(\hat{\beta_{1}}\), se necesitará que:
- Los parámetros verdaderos sean causales.
- MCO deberá proporcionar un buen “estimador” de los valores reales.
Recordemos que \(\hat{\beta_{0}}\) y \(\hat{\beta_{1}}\) son variables aleatorias, no constantes, por lo que un “buen” estimador deberá de ser insesgado y preciso/eficiente.
Pensemos que queremos estimar a través de MCO la siguiente relación lineal entre \(X\) y \(Y\):
\[\begin{equation} Y_{i} = \beta_{0} + \beta_{1}X_{i} + u_{i}, i = 1,...,n \end{equation}\]
MCO será el mejor estimador bajo los siguientes supuestos para todo \(i\):
- $E(u_{i}|X_{i}) = 0 $ para todos los valores de \(X_{i}\).
- \((X_{i}, Y_{i}\) son independientes e identicamente distribuidas.
- Los outliers son poco comunes (curtosis finita).
- \(var(u_{i}|X_{i}) = \sigma_{u}^2\) para todos los valores de \(X_{i}\) (homocedasticidad).
Supuesto 1: \(E(u_{i}|X_{i}) = 0\)
Media condicional igual a cero, se refiere a que la tendencia central de la distribución del término de error \(u_{i}\) no está relacionado a \(X_{i}\), es decir la media de \(u\) es la misma para cualquier valor de \(X\).
Notemos que \(E(u_{i})=0\) no es relevante acá, nos importa \(E(u_{i}|X_{i})\) (esto dado que la media condicional es el mejor predictor lineal en general) sea la misma para todos los valores de \(X\) i.e, \(E(u_{i}) = E(u_{i}|X_{i})\), aunque notemos que \(E(u_{i} = )\) mientras que la regresión tenga una constante
El supuesto de la media condicional igual a cero es más fuerte que el supuesto de no correlación, es decir:
$E(u|X) = 0 corr(u,X) = 0 $ pero
\(corr(u,X) = 0 \notimplies E(u|X) = 0 \newcommand{\notimplies}{\;\not\!\!\!\implies}\)
Notemos que el hecho de que nuestro témino de error \(u\) y \(X\) estén no correlacionados es una condición necesaria pero no suficiente para que el supuesto 1 se cumpla, podemos ilustrar esto de la siguiente forma:
\[\begin{equation} wage = \beta_{0} + \beta_{1}educ + u \end{equation}\]
Necesitamos considerar los factores que serán capturados por \(u . E(u|educ) = o\) implica por ejemplo que \(E(u|trunca\) = \(E(u|grduado\). Esto debido a que existen elemntos que no se capturan con educación en nuestro modelo inicial y serán capturados por \(u\), por lo que \(u\) podría ser la ambición, inteligencia, condiciones de mercado laboral locales, salud, etc. Cualquier cosa que determine \(Y\) y no esté en \(X\) estará en nuestro “término de error”, incluyendo las lurking variables, las cuales serán variables que pueden aparentar tener un nexo causal fuerte de forma errónea.
Si \(X\) está correlacionado con el término de error, el supuesto 1 será violado, y sí \(X\) es asigndo aleatoriamente, entonces no esperaíamos que el término de error esté correlacionado con \(X\), en un RCT el tratamiento es asignado aleatoriamente por lo que el supuesto 1 es plausible.
Por tanto, si el supuesto 1 se cumple entonces esperaríamos que:
Los residuos \(\hat{u}\) siempre estarán no correlacionados a \(X\).
Los residuos \(\hat{u}\) siempre tendrán media cero.
Si corremos la regresión \(\hat{u} = \delta_{0} + \delta_{1}X + \epsilon\) entonces obtenemos $ = = 0 $ siempre.
Los residuos no están correlacionados con \(X\) por construcción.
Nos importa si el error está correacionado o no correlacionado con \(X\).
Por tanto, no podremos testear si \(E(u|X) = 0\) checando si $E(|X) =0 $ dado que esto último siempre es cierto, tampoco podríamos si se cumple \(E(u|X) = 0\) directamente dado que nunca observamos \(u\), por lo que tendremos que decidir si el supuesto 1 es plausible dada la información que sabemos de los datos (por lo que resulta importante pensar en el mecanismo de asignación), otra opción sería observar las relaciones no lineales entre \(\hat{u}\) y \(X\).
Supuesto 2: \((X_{i},Y_{}i)\) son independientes e identicamente distribuidas.
Consideremos que cada unidad de los datos proviene de una realización de la misma población y cada realización es independiente entre sí, notemos que los datos panel generalmente violan esa condición, las series de tiempo también violan ese supuesto regularmente.
Supuesto 3: Los outliers son “poco comunes”.
Es decir \(X\) y \(Y\) tienen curtosis finita, haciendo referencia al cuarto momento central i.e : $ 0 < E(X^4) < $, esta condición también nos sirve para asegurar que tienen una varianza y covarianza finita y bien definida.
Insesgadez de MCO
Consideremos el siguiente modelo especificado para la población:
\[\begin{equation} Y_{i} = \beta_{0} + \beta_{1}X_{i} + u_{i} \end{equation}\]
Si los supuestos 1-3 se cumplen entonces tendremos:
\[\begin{equation} E(\hat{\beta_{0}}) = E(\beta_{0}) \end{equation}\]
\[\begin{equation} E(\hat{\beta_{1}}) = E(\beta) \end{equation}\]
Donde y son estimadores del MCO, lo cual será cierto sin importar la muestra \(n\).
Aproximaciones a la normal (teoría asintótica)
Cuando tenemos una \(n\) lo suficientemente grande, la distribución de \(\hat{\beta}\) se aproxima de forma cercana a la normal, i.e:
\[\begin{equation} \hat{\beta_{1}} N(\beta_{1}, \sigma_{\hat{\beta_{1}}}^2) \end{equation}\]
\[\begin{equation} \sigma_{\hat{\beta_{1}}}^2 = \frac{1}{n} \frac{var [( X - \mu_{X})u ]}{ var(X)^2} \end{equation}\]
Donde \(\mu_{x}\) es la media de la población \(X\).
En este caso usaremos una aplicación del teorema del límite central y nos ayudará a hacer pruebas de hipótesis de forma analítica.
Generalmente se sugiere nuestra muestra \(n >100\) sera suficientemente grande, aunque esto depende en las características de los datos, en algunos casos incluso muestras grandes no serán suficientes dado que no existe suficiente independencia entre los puntos.
Supuesto 4: Homoscedasticidad.
Notemos que aunque es posible \(E(u|X)=0\) pero \(var(u|x)\) puede variar en \(X\), si la varianza de de \(u\) es diferente para diferentes niveles de \(X\), diremos que el término de error es heterocedastico, aunque notemos que si el término de error es homoscedastico entonces nuestros estimadores de MCO serán mejores pero la homoscedasticidad no se requiere para que nuestro MCO sea insesgado
Mejor Estimador Lineal (MEL)
Cuando los supuestos 1-4 se cumplen entonces MCO será el mejor estimador lineal, el cual será 1. el más eficiente/mejor, 2. Lineal (en contraposición con modelos no lineales o de máxima verosimilitud), 3. Insesgado condicionalmente y 4. Estimador.
Variaciones en los estimadores MCO
Recordemos que \(\hat{\beta_{1}}\) es una variable aleatoria, por lo que tendrá una media y si los supuestos 1-3 se cumplen, entonces será insesgado, i.e:
\[\begin{equation} E(\hat{\beta_{1}}) = \beta_{1} \end{equation}\]
Y tiene una desviación estándar de \(\sigma_{\hat{\beta{1}}}\):
\[\begin{equation} \sigma_{\hat{\beta_{1}}}^2 = \frac{1}{n} \frac{var [( X - \mu_{X})u ]}{ var(X)^2} \end{equation}\]
Donde \(n\) es una muestra, el denominador es la varianza de \(X\) al cuadrado.
Sin embargo, es importante preguntarnos: ¿qué determina la precisión de los estimadores MCO?
Pensemos en el caso homoscedastico:
\[\begin{equation} \sigma_{\hat{\beta_{1}}}^2 = \frac{1}{n} \frac{var [( X - \mu_{X})u ]}{ var(X)^2} \end{equation}\]
\[\begin{equation} = \frac{1}{n} \frac{var(X) var(u)]}{ var(X)^2} \end{equation}\]
\[\begin{equation} = \frac{1}{n} \frac{var(u)}{ var(X)^2} = \frac{1}{n} \frac{\sigma_{u}^2}{\sigma_{X}^2} \end{equation}\]
Es decir que tendremos:
\[\begin{equation} = \frac{1}{n} \frac{\sigma_{u}^2}{\sigma_{X}^2} \end{equation}\]
\[\begin{equation} = \frac{1}{\sqrt{n}} \frac{\sigma_{u}}{\sigma_{X}} \end{equation}\]
La precisión de () del estimador MCO es mayor mientras:
- El tamaño de la muestra \(n\) aumente.
- La varianza de \(u\) disminuya.
- La varianza de \(X_{i}\) aumente.
La intuición será que tendremos que observar \(Y\) para un rango más grande de \(X\) para obtener más información.
knitr::include_graphics("Precision_OLS.png")
En el gráfico podemos ver que si bien en ambos casos muestran la verdadera \(\beta\), habrá más varianza en \(X\) del lado derecho, sin embargo ¿cómo podemos saber cuál caso será más preciso? una manera de hacerlo es a través de simulaciones, en tal caso podríamos usar una simulación Monte Carlo como en el siguiente ejemplo:
knitr::include_graphics("MonteCarlo.png")
Podemos ver que al momento de hacer una simulación Monte Carlo con la misma \(\beta\) real, mismo \(n\) y mismo \(\sigma_{u}\) pero con mayor varianza en \(X\) del lado derecho será más preciso dado que obtenemos más información sobre la distribución de \(\hat{\beta_{1}}\)
Medición de ajuste : R cuadrada.
Recordemos que tenemos: \[\begin{equation} Y_{i} = \hat{Y_{i}} + \hat{u_{i}} \end{equation}\]
De lo cual obtenderemos la suma total de los cuadrados: TSS\(= \sum_{i=1}^n(Y_{i} - \bar{Y})^2\), TSS será el valor total de la cantidad de variación en la muestra, notemos que \(var(Y_{i}) = var(\hat{Y_{i}}) + var(\hat{u_{i}})\).
Por otro lado, una vez que tenemos \(\hat{Y_{i}}\) estimado, entonces podemos obtener la suma de cuadrados explicados dada por ESS \(= \sum_{i=1}^n(\hat{Y_{i}} - \bar{Y})^2\) lo cual será la cantidad de variación explicada por \(X\).
Así mismo, tendremos la suma de los residuos cuadrados dado por SSR \(\sum_{i=1}^n \hat{U}^2\) que será la cantidad de variación no explicada por \(X\), notemos que TSS = ESS+SSR.
A partir de lo anterior podemos definir la R cuadrada \(R^2\) la cual se define como la fracción de la varianza muestral en \(Y_{i}\) que está explicado por \(X_{i}\), es decir:
\[\begin{equation} R^2 = \frac{ESS}{TSS} \end{equation}\]
\[\begin{equation} R^2 = 1 - \frac{SSR}{TSS} \end{equation}\]
Donde la \(R^2\) oscila entre 0 y 1, por lo que:
- \(R^2 = 1 \iff SSR=0\)
- \(R^2 = 0 \iff SSR=0\)
Por lo que una cantidad más cercana a 1 un mejor ajuste de nuestro modelo o bien una proporción más grande de la variación explicada.
Medición de ajuste : Error estándar de la regresión
El error estándar de la regresión es una medida del error de predicción y usado para estimar la varianza de los coeficientes, un estimador de la desviación estándar de \(u_{i}\) es el error estándar de la regresión (SER) dado como:
\[\begin{equation} SER = \frac{1}{n-2} \sum_{i=1}^n \hat{u_{i}^2} = \frac{1}{n-2} SSR \end{equation}\]
El SER mide la típica desviación entre \(Y_{i}\) y \(\hat{Y_{i}}\) y \(n-2\) será los níveles de libertad (tamaño de la muestra menos el número de parámetros estimados).
knitr::include_graphics("Fitted.png")
Notemos que los datos del lado
derecho son los mismos que los del lado izquierdo con el error reducido
(con varianza menor de \(u\)), por lo
que la gráfica del lado derecho tendrá una mayor \(R^2\) dado que una fracción de la variación
en \(Y\) es explicada por \(X\).
MCO : Inferencia - Hipótesis Alternativa/ Hipótesis Nula
Vamos a concentrarnos en dos conceptos importantes:
- \(H_0\) : Hipoótesis Nula refiriéndose a que el efecto no es \(estádisticamente \ significativo\) o bien que no hay efecto.
- \(H_A\) : Hipótesis Alternativa, refiriéndose a que el efecto es \(estádisticamente \ significativo\).
Generalmente buscamos desacreditar la hipótesis nula como objetivo de la investigación y hacer inferencia estadística sobre hipótesis relativo a los parámetros de la población \(\beta\), de tal forma podemos decir que:
\(H_0: \beta_{1} = \theta\) vs \(H_{A}: \beta_{1} \neq \theta\) (two-sided)
\(H_0: \beta_{1} = \theta\) vs \(H_{A}: \beta_{1} < \theta\) (one-sided)
Generalmente, se usan los test two-sided de forma más apropiada.
Intervalos de confianza
Imaginemos que tenemos \(\hat{\beta} = 0.084\) y el error estándar de \(\hat{\beta}\) (SE(\(\beta\))) = 0.014.
Con la información anterior podríamos construir elintervalo de confianza a 95%a para \(\beta\) (usando una aproximación normal) como: \([\hat{\beta}-1.96 *SE(\hat{\beta}), \hat{\beta}+1.96*SE(\hat{\beta})] = [0.0567, 0.114]\) la primera parte nos dará la cola final y la segunda parte nos dará la cola inicial.
Lo anterior se puede traducir también como que el 95% de las muestras posibles, el Intervalo de Confianza (CI) incluirá el valor verdadero.
Para cualquier \(\theta\) dentro del Intervalo de Confianza (CI) se incluirá el valor verdadero, para cualquier \(\theta\) dentro del CI no se rechazará $H_{0} : $ a un nivel del 5%.
Estadístico t
Será ampliamente utilizado, el test default en Stata será contra cero (no efecto) efecto nulo, el estadístico t de $H_{0}: _{1} = 0 $, tomando el ejemplo anterior: $ = 6.0 $
Pero podemos hacer el mismo test contra cualquier otra hipótesis nula, el estadístico t de $H_{0}: _{1} = 0.1 $, tomando el ejemplo anterior: $ = -1.14 $
Podemos recordar que si W es una variable normal, si le restamos la media y dividimos sobre la desviación estándar nos lleva a una variable normal estándar, es decir:
\[\begin{equation} W - N(\mu_{W}, \sigma_{W}^2) \rightarrow \frac{W - \mu_{W}}{\sigma_{W}} - N(0,1) \end{equation}\]
Si la desviación estándar es estimada en lugar de estar dada, entonces al momento de estandariar nos llevará a una distribución t
\[\begin{equation} W - N(\mu_{W}, \sigma_{W}^2) \rightarrow \frac{W - \mu_{W}}{\hat{\sigma_{W}}} - t_{n-k-1} \end{equation}\]
Y cuando n-k-1 es grande, entonces \(t\) es aproximadamente normal.
knitr::include_graphics("Est_t.png")
knitr::include_graphics("Est_t2.png")
knitr::include_graphics("Est_t3.png")
knitr::include_graphics("Est_t4.png")
Podemos interpretar los grados de libertad, pero ¿cuándo será n-k-1 suficientemente grande? podemos ver en el gráfico anterior cuándo la distribución t empata con la distribución normal, mientras que cuando aumentamos los grados de libertad de 1 a 5, 10 y 30 se acerca más a la distribución normal.
Así como el MCO será nuestro estimador estándar del regresor, el \(t-stat\) será nuestro test estadístico estándar.
Sobre la hipótesis nula de que \(\beta\) es igual a alguna constante \(\theta\) (\(H_{0}: \beta = \theta\)) podemos formular el estadístico t:
\[\begin{equation} \frac{\hat{\beta - \theta}}{\hat{\sigma_{\hat{\beta}}}} -.- t_{n-k-1} \end{equation}\]
Si el estadístico t es muy pequeño, no podremos rechazar la hipótesis nula.
También existe una “regla del pulgar” la cual implica que \(|t-stat| >2\) implica que podemos rechazar la hipósteis nula.
knitr::include_graphics("Est_t5.png")
No rechazamos la hipótesis nula \(H_{o}\) si \(|t-stat|<2\) dado que el valor absoluto de ese estadístico t es menor que 2.
knitr::include_graphics("Est_t6.png")
En este caso sí rechazamos la hipótesis nula \(H_{o}\) si \(|t-stat|>2\) dado que el valor absoluto de ese estadístico t es mayor que 2.
Los resultados anteriores se pueden observar en la siguiente tabla, misma que generalmente se muestra después de correr un modelo estadístico:
knitr::include_graphics("Tabla_est.png")
Podemos observar como la constante será 70.3973 lo cual significa que en esta muestra en particular tendremos individuos que son aproximadamente 70.39 inch de altura cuando female (mujer) = 0, ¿qué pasa cuando female =1? utilizaremos el valor de -10.322 menos cero, la hipótesis nula en ese sentido sería que no existe una diferencia entre las alturas de hombres y mujeres en la muestra, y será dividio sobre el error estándar (SE) dado por 1.699
Stata automaticamente reporta \(H_{0}: \beta_{1} = 0\), en el caso que hemos estado discutiendo sería:
\[\begin{equation} \frac{-10.33 -0}{1.699} = -6.08 \end{equation}\]
Recordemos que el test nos muestra si el estimador del coeficiente es “estadísticamente significativo diferente de cero”, dado que tomaremos el valor absoluto, por lo que que sea negativo no impacta realmente, notemos que es mayor a 2 por lo que podemos refutar fácilmente la hipótesis nula en este caso, este test también sugiere que el coeficiente estimado es estadísticamente significativo, dado que será diferente de 0 donde 0 es nuestra hipótesis nula, generalmente se toman estos resultados como muestras de patrones descriptivos que pueden sugerir causalidad.
Valores \(p\)
Los valores \(p\) se definen como la probabilidad de obtener una \(\hat{\beta}\) que sea lo suficientemente lejana al valor nulo como del resultado observado, si el nulo es verdadero entonces definimos:
\[\begin{equation} p - value \aprox 2\Phi(-|t|) \end{equation}\]
El resultado anterior usa una aproximación de la normal estándar para la distribución t con diferentes niveles de libertad.
Parecido al estadístico t, lo anterior depende del valor nulo en cuestión:
\(p\)-value (para \((H_{0}: \beta = 0) \aprox 2\Phi (-6.08) < 0.00001\))
\(p\)-value (para \((H_{0}: \beta = -9) \aprox 2\Phi (0.782) = 0.439\))
MCO: Causalidad
Si bien es algo común encontrar correlacoones entre diferentes variables, no quiere decir que tods esas correlaciones nos lleven a causalidad, supongamos que hay una relación lineal causal: \(X\) causa a \(Y\).
Si estimamos por MCO $ X = {0} + {1}Y$ podemos generar una \(\beta_{1}\) diferente de cero y tiene una interpretación inversa.
Incluso si tuviésmos la dirección de la causadlidad de forma correcta, MCO podría danos un estimador sesgado para \(\beta_{1}\) en la mayoría de los casos.
Bajo circunstancias especiales MCO nos puede dar un etimador “bueno” de la relación causal entre \(X\) y \(Y\).
Cuando los datos se generan de un RCT realizado correctamente, al momento de aplicar MCO obtendremos un estimador insesgado para una relación causal dado que de inicio tenemos correctamente identificado nuestro nexo causal y también se satisfacen las condiciones del Teorema de Gauss-Markov.
Podemos pensar en el caso en el que tenemos un tratamiento bivariado, los valores predecidos (esperanas condicionadas estimada) si \(X\) es una variable binaria se definen de la siguiente forma:
Definamos una variable dummy como =1 si el enunciado es verdadero, 0 en otro caso. Por ejemlo, una variable dummy de tratamiento será 1 si \(i\) es tratado y o en otro caso (si \(i\) forma parte del grupo de control).
\[\begin{equation} \hat{E}(Y_{i}|X_{i} = 1) = \hat{\beta_{0}} + \hat{\beta_{1}} x 1 = \hat{\beta_{0}} + \hat{\beta_{1}} \end{equation}\]
\[\begin{equation} \hat{E}(Y_{i}|X_{i} = 0) = \hat{\beta_{0}} + \hat{\beta_{1}} x 0 = \hat{\beta_{0}} \end{equation}\]
library(ggplot2)
library(NHANES)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorscolnames(NHANES)## [1] "ID" "SurveyYr" "Gender" "Age"
## [5] "AgeDecade" "AgeMonths" "Race1" "Race3"
## [9] "Education" "MaritalStatus" "HHIncome" "HHIncomeMid"
## [13] "Poverty" "HomeRooms" "HomeOwn" "Work"
## [17] "Weight" "Length" "HeadCirc" "Height"
## [21] "BMI" "BMICatUnder20yrs" "BMI_WHO" "Pulse"
## [25] "BPSysAve" "BPDiaAve" "BPSys1" "BPDia1"
## [29] "BPSys2" "BPDia2" "BPSys3" "BPDia3"
## [33] "Testosterone" "DirectChol" "TotChol" "UrineVol1"
## [37] "UrineFlow1" "UrineVol2" "UrineFlow2" "Diabetes"
## [41] "DiabetesAge" "HealthGen" "DaysPhysHlthBad" "DaysMentHlthBad"
## [45] "LittleInterest" "Depressed" "nPregnancies" "nBabies"
## [49] "Age1stBaby" "SleepHrsNight" "SleepTrouble" "PhysActive"
## [53] "PhysActiveDays" "TVHrsDay" "CompHrsDay" "TVHrsDayChild"
## [57] "CompHrsDayChild" "Alcohol12PlusYr" "AlcoholDay" "AlcoholYear"
## [61] "SmokeNow" "Smoke100" "Smoke100n" "SmokeAge"
## [65] "Marijuana" "AgeFirstMarij" "RegularMarij" "AgeRegMarij"
## [69] "HardDrugs" "SexEver" "SexAge" "SexNumPartnLife"
## [73] "SexNumPartYear" "SameSex" "SexOrientation" "PregnantNow"# Create bar plot for Home Ownership by Gender
ggplot(NHANES, aes(x = Gender, fill = HomeOwn)) +
# Set the position to fill
geom_bar(position = "fill") +
ylab("Relative frequencies")
# Density plot of SleepHrsNight colored by SleepTrouble
ggplot(NHANES, aes(x = SleepHrsNight, color = SleepTrouble)) +
# Adjust by 2
geom_density(adjust = 2) +
# Facet by HealthGen
facet_wrap(~ HealthGen)## Warning: Removed 2245 rows containing non-finite outside the scale range
## (`stat_density()`).
#Using the NHANES dataset, let's investigate the relationship between gender and home ownership. Remember, more information #about the dataset can be found here: NHANES.
# Subset the data: homes
homes <- NHANES %>%
# Select Gender and HomeOwn
select(Gender, HomeOwn) %>%
# Filter for HomeOwn equal to "Own" or "Rent"
filter(HomeOwn %in% c("Own", "Rent"))Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.